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Abstract This is a review of the theoretical background, experimental techniques, and phe- 
nomenology of what is called the "Glauber Model" in relativistic heavy ion physics. This model 
is used to calculate "geometric" quantities, which are typically expressed as impact parameter 
(&), number of participating nucleons (N pax t) and number of binary nucleon-nucleon collisions 
(N co u). A brief history of the original Glauber model is presented, with emphasis on its devel- 
opment into the purely classical, geometric picture that is used for present-day data analyses. 
Distinctions are made between the "optical limit" and Monte Carlo approaches, which are often 
used interchangably but have some essential differences in particular contexts. The methods 
used by the four RHIC experiments are compared and contrasted, although the end results 
are reassuringly similar for the various geometric observables. Finally, several important RHIC 
measurements are highlighted that rely on geometric quantities, estimated from Glauber cal- 
culations, to draw insight from experimental observables. The status and future of Glauber 
modeling in the next generation of heavy ion physics studies is briefly discussed. 
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1 Introduction 

Ultra-relativistic collisions of nuclei produce the highest multiplicities of outgoing 
particles of all subatomic systems known in the laboratory. Thousands of particles 
are created when two nuclei collide head-on, generating event images of dramatic 
complexity compared with proton-proton collisions. The latter is a natural point 
of comparison, as nuclei are themselves made up of nucleons, i.e. protons and 
neutrons. Thus, it is natural to ask just how many nucleons are involved in any 
particular collision or, more reasonably, in a sample of selected collisions. It is 
also an interesting to consider other ways to characterize the overlapping nuclei, 
e.g. their shape. 

While this problem would seem intractable, with the femtoscopic length scales 
involved precluding direct observation of the impact parameter (6) or number of 
participating nucleons (-/V par t) or binary nucleon-nucleon collisions (N co \\), the- 
oretical techniques have been developed to allow estimation of these quantities 
from experimental data. These techniques, which consider the multiple- scatter- 
ing of nucleons in nuclear targets, are generally referred to as "Glauber Models" 
after Roy Glauber. Glauber pioneered the use of quantum mechanical scattering 
theory for composite systems, describing non trivial effects discovered in cross 
sections for proton-nucleus (p+A) and nucleus- nucleus (A+B) collisions at low 
energies. 

Over the years, a variety of methods were developed for calculating geomet- 
ric quantities relevant for p+A and A+B collisions dUElE]). Moreover, a wide 
variety of experimental methods were tried to connect distributions of measured 
quantities to the distributions of geometric quantities. This review is an attempt 
to explain how the RHIC experiments grappled with this problem, and largely 
succeeded, both because of good-sense experimental and theoretical thinking. 

Heavy ion physics entered a new era with the turn-on of the RHIC collider in 
2000 (4). Previous generations of heavy ion experiments searching for the Quark 
Gluon Plasma (QGP) had focused on particular signatures suggested by theory. 
From this strategy, experiments generally measured observables in different re- 
gions of phase space, or focused on particular particle types. RHIC experiments 
were designed in a comprehensive way, with certain regions of phase space being 
covered by multiple experiments. This allowed systematic cross checks of vari- 
ous measurements between experiments, increasing the credibility of the separate 
results dUElElED. 

Among the most fundamental observables shared between the experiments were 
those relating to the geometry of the collision. Identical Zero-Degree Calorime- 
ters ([9]) were installed in all experiments to estimate the number of spectator 
nucleons, and all experiments had coverage for multiplicities and energy mea- 
surements over a wide angular range. This allowed a set of systematic studies of 
centrality in d+Au and A+A collisions, providing one of the first truly extensive 
datasets all of which characterized by geometric quantities. 

This review will cover the basic information a newcomer to the field should 
know to understand how centrality is estimated by the experiments, and how 
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this is related to nucleus- nucleus collisions. Section [2] will discuss the history of 
the Glauber model by reference to the theoretical description of nucleus-nucleus 
collisions. Section [3] will discuss how experiments measure centrality by a variety 
of methods and relate that by a simple procedure to geometrical quantities. Sec- 
tion H] will illustrate the relevance of various geometrical quantities by reference 
to actual RHIC data. These examples will illustrate how a precise quantita- 
tive grasp of the geometry allows qualitatively new insights into the dynamics 
of nucleus-nucleus collisions. Finally, section [5] will assess the current state of 
knowledge and point to new directions in our understanding of nuclear geometry 
and its impact on present and future measurements. 
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2 Theoretical foundations of Glauber modeling 
2.1 A brief history of the Glauber model 

The Glauber model was developed to address the problem of high energy scat- 
tering with composite particles. This was of great interest to both nuclear and 
particle physicists, who have both benefited from the insights of Glauber in the 
1950's. In his 1958 lectures, Glauber presented his first collection of various 
papers and unpublished work from the 1950's ([2]). Some of this was updated 
work started by Moliere in the 1940's, but much of it was in direct response to 
the new experiments involving protons, deuterons and light nuclei. Up to that 
point, there had been no systematic calculations treating the many-body nuclear 
system either as a projectile or target. Glauber's work put the quantum the- 
ory of collisions of composite objects on a firm basis, and provided a consistent 
description of experimental data for protons incident on deuterons and larger 
nuclei (fTHfTUl) . Most striking were the observed dips in the elastic peaks whose 
position and magnitude were predicted by Glauber's theory by Czyz and Lesniak 
in 1967 (fT2|) . 

It was only in the 1970's when high energy beams of hadrons and nuclei were 
systematically scattered off of nuclear targets. Glauber's work was found to have 
utility in describing total cross sections, for which "factorization" relationships 
were found (e.g. o\ B ~ (Jaa&bb) (|13pi4p . Maximon and Czyz applied the theory 
in its most complete form to p+A and A+B collisions in 1969, focusing mainly 
on elastic collisions (JT). Finally, Bialas, Bleszynski, and Czyz (115l l3|) applied 
Glauber's approach to inelastic nuclear collisions, after they had already applied 
their "wounded nucleon model" to hadron-nucleus collisions This is essentially 
the bare framework of the traditional "Glauber Model" , with all of the quantum 
mechanics reduced to its simplest form f)16|) . The main remaining feature of the 
original Glauber calculations is the "optical limit" , used to make the full multiple 
scattering integral numerically tractable. 

The approach of Bialas et al. (15) introduced the basic functions used in more 
modern language, including the thickness function and a prototype of the nuclear 
overlap function Tab- This paper also introduced the convention of using the 
optical limit for analytic and numerical calculations, despite full knowledge that 
the "real" Glauber calculation is an 2{A + B + l)-dimensional integral over the 
impact parameter and each of the nucleon positions. A similar convention exists 
in most theoretical calculations of geometrical parameters to this day. 

With the advent of desktop computers, the "Glauber Monte Carlo" (GMC) 
approach emerged as a natural framework for use by more realistic particle pro- 
duction codes (I17 pi8h . The idea was to model the nucleus in the simplest way, as 
uncorrelated nucleons sampled from measured density distributions. Two nuclei 
could be arranged with a random impact parameter b and projected onto the 
x-y plane. Then interaction probabilities could be applied by using the relative 
distance between two nucleon centroids as a stand-in for the measured nucleon- 
nucleon inelastic cross section. 

The GMC approach was first applied to high energy heavy ion collisions in the 
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Figure 1: Density distributions for nuclei used at RHIC (a) and distribution 
of the proton-neutron distance in the deuteron as given by the Hulthen wave 
function (b). 



HIJET model (|1T|) and has shown up in practically all A+A simulation codes. 
This includes HIJING (32]), VENUS (20), RQMD (JHj), and all models which 
require specific production points for the different sub-processes in a nucleus- 
nucleus collision, rather than just aggregate values. 



2.2 Inputs to Glauber Calculations 

In all calculations of geometric parameters using a Glauber approach, some ex- 
perimental data must be given as model inputs. The two most important are the 
nuclear charge densities measured in low-energy electron scattering experiments 
and the energy dependence of the inelastic nucleon-nucleon cross section. 
2.2.1 Nuclear Charge Densities The nucleon density is usually param- 
eterized by a Fermi distribution with three parameters: 

l + w(r/R) 2 

p( r ) = Po 7— I 1 ) 

1 + exp (r=S) 

where po corresponds to the nucleon density in the center of the nucleus, R 
corresponds to the nuclear radius, a to the "skin depth" and w characterizes 
deviations from a spherical shape. For 197 Au (R = 6.38 fm; a = 0.535 fm; w = 0) 
and 63 Cu (i? = 4.20641 fm; a = 0.5977 fm; w = 0), the nuclei so far employed at 
RHIC, p{r ) I pq is shown in Fig. [T^, with the Fermi distribution parameters as given 
in ref. (22,23). In the Monte Carlo procedure the radius of a nucleon is randomly 
drawn from the distribution 47rr 2 p(r) (where the absolute normalization is of 
course irrelevant). 
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Figure 2: Inelastic nucleon- nucleoli cross section as parameterized by Pythia 
(|34ll35|) (solid line) along with data on total and elastic nucleon-nucleon cross 
sections as function of y/s (36). The stars indicate nucleon-nucleon cross section 
used for Glauber Monte Carlo calculations at RHIC (c^ = 32.3, 35.6, 40, and 
42 mb at y/s^ = 19.6, 62.4, 130, and 200 GeV) 



At RHIC, effects of cold nuclear matter have been studied with the aid of 
d+Au collisions. In the Monte Carlo calculations the deuteron wave function 
was represented by the Hulthen form (J2] 



1 y/ab(a + b) e~ ar <™ - e~ br ^ 



2ir b — a r pn 

with parameters a = 0.228 fm" 1 and b = 1.18 fm -1 (26). The variable r pn in 
Eq. [2] denotes the distance between the proton and the neutron. Accordingly, 
r pn was drawn from the distribution p(r pn ) = ^irr^cf) 2 (r pn ) which is shown in 
Fig. Hb. 

2.2.2 Inelastic Nucleon-Nucleon Cross Section (ct^i) In the con- 
text of high energy nuclear collisions, we are typically interested in multiparticle 
nucleon-nucleon processes. As the cross section involves processes with low mo- 
mentum transfer, it is impossible to calculate this using perturbative QCD. Thus, 
the measured inelastic nucleon-nucleon cross section (<r?^j) is used as an input, 
and provides the only non-trivial beam-energy dependence for Glauber calcula- 
tions. From y/s^ ~ 20 GeV (CERN SPS) to y/a^ = 200 GeV (RHIC) <t?£J 
increases from ~ 32 mb to ~ 42 mb as shown in Fig. [2j Diffractive and elastic 
processes, which are typically ignored in high energy multiparticle nuclear colli- 
sions, are generally active out to large impact parameters, and thus require full 
quantum mechanical wave functions. 
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Figure 3: Schematic representation of the Optical Glauber Model geometry, with 
transverse (a) and longitudinal (b) views. 

2.3 Optical-limit Approximation 

The Glauber Model views the collision of two nuclei in terms of the individual 
interactions of the constituent nucleons (see, e.g., Ref. (|27p ). In the optical limit, 
the overall phase shift of the incoming wave is taken as a sum over all possible 
two-nucleon (complex) phase shifts, with the imaginary part of the phase shifts 
related to the nucleon-nucleon scattering cross section through the optical theo- 
rem ([28j[29]). The model assumes that at sufficiently high energies, these nucleons 
will carry sufficient momentum that they will be essentially undeflected as the 
nuclei pass through each other. It is also assumed that the nucleons move inde- 
pendently in the nucleus and that the size of the nucleus is large compared to the 
extent of the nucleon-nucleon force. The hypothesis of independent linear tra- 
jectories of the constituent nucleons makes it possible to develop simple analytic 
expressions for the nucleus-nucleus interaction cross section and for the number 
of interacting nucleons and the number of nucleon-nucleon collisions in terms of 
the basic nucleon-nucleon cross section. 

Consider Fig. [3l Two heavy-ions, "target" A and "projectile" B are shown 
colliding at relativistic speeds with impact parameter b (for colliding beam ex- 
periments the distinction between the target and projectile nuclei is a matter of 
convenience). We focus on the two flux tubes located at a displacement s with 
respect to the center of the target nucleus and a distance s — b from the center 
of the projectile. During the collision these tubes overlap. The probability per 
unit transverse area of a given nucleon being located in the target flux tube is 
Ta (s) = / Pa(&, ZA)dzA, where pa (s, za) is the probability per unit volume, nor- 
malized to unity, for finding the nucleon at location (s, za)- A similar expression 
follows for the projectile nucleon. The product Ta (s) Tb (s — b) d 2 s then gives 
the joint probability per unit area of nucleons being located in the respective 
overlapping target and projectile flux tubes of differential area d 2 s. Integrating 
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this product over all values of s defines the "thickness function" T (b), with 

f AB (b) = jf A (s) f B (s - b)d 2 s. (3) 



Notice that T (b) has the unit of inverse area. We can interpret this as the 
effective overlap area for which a specific nucleon in A can interact with a given 
nucleon in B. The probability of an interaction occurring is then T (b) cxj^, where 
°inei ^ s ^ ne nucleon-nucleon inelastic cross section. Elastic processes lead to very 
little energy loss and are consequently not considered in the Glauber-model cal- 
culations. Once the probably of a given nucleon-nucleon interaction has been 
found, the probably of having n such interactions between nucleus A (with A 
nucleons) and B (with B nucleons) is given as a binomial distribution, 



P(n,b) 



I AB 

V n 



f AB (b) u. 



NN 
incl 



1 - f AB (b) a. 



NN 
incl 



AB-n 



(4) 



where the first term is the number of combinations for finding n collisions out 
of AB possible nucleon-nucleon interactions, the second term the probability for 
having exactly n collisions, and the last term is the probability of exactly AB — n 
misses. 

Based on this probability distribution, a number of useful reactions quantities 
can be found. The total probability of an interaction between A and B is 



'2A.+B 



incl 



db 2 



— Pincl 



AB 

E 

71=1 



P (n, b) 



1 



1 - Tab (b) 



AB 



(5) 



The vector impact parameter can be replaced by a scalar distance if the nuclei 
are not polarized. In this case, the total cross section can be found as 



A+B 
inel 



2Mb { 1 



1 - Tab (b) a 



NN 
inel 



AB 



(6) 



The total number of nucleon-nucleon collisions is 



iVcoii (6) 



AB 

E 

71=1 



nP (n, b) = ABTab (6) ° 



NN 
incl 



(7) 



using the result for the mean of a binomial distribution. The number of nucleons 
in the target and projectile nuclei that interact is called either the "number of 
participants" or the "number of wounded nucleons" . The number of participants 
(or wounded nucleons) at impact parameter b is given by (|15 pop 

si 



iV part (b) = A 



B 



1 



T A (s 

T B (s - b) { 1 



1 - T B (s 



K^i ^NN 
b ) ^incl 



1 - T A (s) aZ 



d 2 s+ 



d 2 s, 



(8) 



where it can be noted that the integral over the bracketed terms give the respec- 
tive inelastic cross sections for nucleon-nucleus collisions: 

A(B) 



A(B) 
'inel 



d 2 s{ 1 



1 - a fr^\TA{B) (s) 



(9) 
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Figure 4: Glauber Monte Carlo event (Au+Au at ^/snn = 200 GeV with impact 
parameter 6 = 6 fm) viewed in the transverse plane (left panel) and along the 
beam axis (right panel). The nucleons are drawn with a radius \J a^ l /ir/2. 
Darker disks represent participating nucleons. 



The optical form of the Glauber theory is based on continuous nucleon density 
distributions. The theory does not locate nucleons at specific spatial coordinates, 
as is the case for the Monte Carlo formulation that is discussed in the next section. 
This difference between the optical and Monte Carlo approaches can lead to subtle 
differences in calculated results, as will be discussed below. 

2.4 Glauber Monte Carlo approach 

The virtue of the Monte Carlo approach for the calculation of geometry related 
quantities like (iV part ) and (iV co ii) is its simplicity. Moreover, it is possible to 
simulate experimentally observable quantities like the charged particle multi- 
plicity and to apply similar centrality cuts as in the analysis of real data. In 
the Monte Carlo ansatz the two colliding nuclei are assembled in the computer 
by distributing the A nucleons of nucleus A and B nucleons of nucleons B in 
three-dimensional coordinate system according to the respective nuclear density 
distribution. A random impact parameter b is then drawn from the distribution 
dcx/d6 = 2-Kb. A nucleus-nucleus collision is treated as a sequence of indepen- 
dent binary nucleon-nucleon collisions, i.e., the nucleons travel on straight-line 
trajectories and the inelastic nucleon-nucleon cross-section is assumed to be inde- 
pendent of the number of collisions a nucleon underwent before. In the simplest 
version of the Monte Carlo approach a nucleon-nucleon collision takes place if 
their distance d in the plane orthogonal to the beam axis satisfies 

^<V^> (io) 

where crP^j is the total inelastic nucleon-nucleon cross-section. As an alterna- 
tive to the black-disk nucleon-nucleon overlap function, e.g., a Gaussian overlap 
function can be used ([3T|) . 
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Figure 5: (left) Total cross section, calculated in the optical approximation and 
with a Glauber Monte Carlo, both with identical nuclear parameters, as a function 
of 

^ineL the nucleon-nucleon inelastic cross section, (right) N co \\ and Ap art as a 
function of impact parameter, calculated in the optical approximation (lines) and 
with a Glauber Monte Carlo (symbols). The two are essentially identical out to 
b = 2R A . 



An illustration of a Glauber Monte Carlo event for a Au+Au collision with 
impact parameter b = 6 fm is shown in Fig. SJ The average number of partic- 
ipating nucleons and binary nucleon-nucleon collisions and other quantities are 
then determined by simulating many nucleus-nucleus collisions. 

2.5 Differences between Optical and Monte Carlo Approaches 

It is not always remembered that the various integrals used to calculated physical 
observables in the "Glauber Model" are predicated on a particular approximation 
called the optical limit. This limit assumes that scattering amplitudes can be 
described by an eikonal approach, where the incoming nucleons see the target as 
a smooth density. This approach captures many features of the collision process, 
but does not completely capture the physics of the total cross section. Thus, 
it tends to lead to distortions in the estimation of A r part and N co \\ compared to 
similar estimations using the Glauber Monte Carlo approach. 

This can be seen by simply looking at the relevant integrals. The full (2A+2B+2)- 
dimensional integral to calculate the total cross section is f| 1 5 1> : 

o AB = J d 2 bj d 2 sf ■ ■ ■ d 2 sid 2 sf ■ ■ ■ d 2 s% x 

f A {s£) ■ ■ ■ f A (si)f B (sf) ■ ■ ■ f B (sl) x (11) 

ji-nn[i-*(b-^+sf)]j 

I 3=1 i=l J 
where <r(s) is normalized such that Jd 2 sa(s) = cr-^, while the optical limit 
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Figure 6: a) Dependence of N pait and iV co ii on cr^ for central (b < 2 fm) and 
minimum bias (MB) events for the 197 Au + 197 Au reaction, b) Dependence of 
these values on the system size, with = 42 mb. 



version of the same calculation is (cf. Equation [6]) : 

oab = Jd%{l-[l- al*f AB (b)] AB } (12) 

These expressions are generally expected to be the same for large A (and B) 
and/or when is sufficiently small (115h . The main difference between the 
two is that many terms in the full calculation are missing in the optical limit 
calculation. These are the terms which describe local density fluctuations event- 
by-event. Thus, in the optical limit, each nucleon in the projectile sees the 
oncoming target ciS 3j smooth density. 

One can test this interpretation to first order by plotting the total cross section 
as a function of crP^j for an optical limit calculation as well as a GMC calculation 
with the same parameters, as shown in the left panel of Fig. [5j One sees that as 
the nucleon-nucleon cross section becomes more point-like, the optical and GMC 
cross sections converge. This confirms the general suspicion that it is the ability 
of GMC calculations to introduce "shadowing" corrections that reduces the cross 
section relative to the optical calculations (jlip 

And yet, when calculating simple geometric quantities like N paTt and iV co ii as 
a function of impact parameter, there is little difference between the two calcula- 
tions, as shown in the right panel of Fig. [5j The only substantial difference comes 
at the highest impact parameters, something which will be discussed in Section 
13.4.11 Fluctuations are also sensitive to this difference, but there is insufficient 
room in this review to discuss more recent developments (|32p , 

2.6 Glauber Model Systematics 

As discussed above, the Glauber model depends on the nucleon-nucleon inelastic 
cross <7jj[5 section and the geometry of the interacting nuclei. In turn, cr^S 
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Figure 7: Geometric scaling behavior of N co \\ as discussed in the text. The 
calculations are done with = 42 mb. 



depends on the energy of the reaction, as shown in Fig. (2J and the geometry 
depends on the number of nucleons in the target and projectile. 

Figure [6^ shows the effect of changing cr-^ on the calculated values of the num- 
ber of participant (-/Vp art ) and the number of collisions (N co \\) for 197 Au + 197 Au 
reaction over the range of values relevant at RHIC. The secondary axis shows 
the values of - v /5nn corresponding to the <r^ values. The values are shown for 
central events, with impact parameter b < 2 fm, and for a minimum bias throw 
of events. The error bars, which only extend beyond the symbol size for the 
Acoii(fr < 2 fm) results, are based on an assumed uncertainty of 0";^] at a given 
energy of 3 mb. In general, one finds that the Glauber calculations show only a 
weak energy dependence over the energy range covered by the RHIC accelerator. 

Fig. [6b shows dependence of N part and A^ou on the system size for central and 
MB events, with values calculated for identical particle collisions of the indicated 
systems at a fixed value of tr?^ = 42 mb (corresponding to y / «NN = 200 GeV). 
Since the Glauber Model is largely dependent on the geometry of the colliding 
nuclei, some simple scalings can be expected for Ap art and N co \\. 

Ap ar t should scale with the volume of the interaction region. In Fig. [6] this is 
seen by the linear dependence of AT part for central collisions on A, where the com- 
mon volume V of the largely overlapping nuclei in central collisions is proportional 
to A for a saturated nuclear density. In a collisions of two equal nuclei (A+A) 

1/3 

the average number of collisions per participant scales as the length l z oc A^ part of 
the interaction volume along the beam direction so that the number of collisions 
roughly follows 

iV coll oc iV^ (13) 

independent of the size of the nuclei. This scaling relationship is demonstrated 
in Fig. [7] where N co \\/N^ t is plotted as a function of A^ part for the systems shown 
in Fig. [6)3. The result is confirmed experimentally in ref. (|33p where the Cu + Cu 
and Au + Au systems are compared. The geometric nature of the Glauber Model 
is evident. 
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3 Relating the Glauber Model to Experimental Data 

Unfortunately, neither iV par t nor N co \\ can be directly measured in a RHIC exper- 
iment. Mean values of such quantities can be extracted for classes of (iV ev t) mea- 
sured events via a mapping procedure. Typically a measured distribution (e.g., 
diVevt/diVch) is mapped to the corresponding distribution obtained from phe- 
nomenological Glauber calculations. This is done by defining "centrality classes" 
in both the measured and calculated distributions and then connecting the mean 
values from the same centrality class in the two distributions. The specifics of this 
mapping procedure differ both between experiments as well as between collision 
systems within a given experiment. Herein we briefly summarize the principles 
and various implementations of centrality definition. 



3.1 Methodology 




Figure 8: A cartoon example of the correlation of the final state observable 
JVch with Glauber calculated quantities (b, N part ). The plotted distribution and 
various values are illustrative and not actual measurements (T. Ullrich, private 
communication) . 



The basic assumption underlying centrality classes is that the impact param- 
eter b is monotonically related to particle multiplicity, both at mid and forward 
rapidity. For large b events ("peripheral") we expect low multiplicity at mid- 
rapidity, and a large number of spectator nucleons at beam rapidity, whereas 
for small b events ("central") we expect large multiplicity at mid-rapidity and a 
small number of spectator nucleons at beam rapidity (Figure [8]) . In the simplest 
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Figure 9: Pseudorapidity coverage of the centrality selection detectors of the 4 
RHIC experiments. Top to bottom: BRAHMS, PHENIX, PHOBOS, STAR. 



case, one measures the per-event charged particle multiplicity (diV evt /diV c h) for 
an ensemble of events. Once the total integral of the distribution is known, cen- 
trality classes are defined by binning the distribution based upon the fraction 
of the total integral. The dashed vertical lines in Figure [5] represent a typical 
binning. The same procedure is then applied to a calculated distribution, often 
derived from a large number of Monte Carlo trials. For each centrality class, the 
mean value of Glauber quantities (e.g., {N par t}) for the Monte Carlo events within 
the bin (e.g., 5-10%) is calculated. Potential complications to this straightfor- 
ward procedure arise from various sources: event selection, uncertainty in the 
total measured cross section, fluctuations in both the measured and calculated 
distributions, and finite kinematic acceptance. 

3.1.1 Event Selection All four RHIC experiments share a common de- 
tector to select minimum bias (MB) heavy ion events. The Zero Degree Calorime- 
ters (ZDCs) are small acceptance hadronic calorimeters with angular coverage of 
8 < 2 mrad with respect to the beam axis ((9|). Situated behind the charged par- 
ticle steering DX magnets of RHIC, the ZDCs are primarily sensitive to neutral 
spectators. For Au+Au collisions at ^/snn = 130 GeV and above, the ZDCs are 
~ 100% efficient for inelastic collisions, thus providing an excellent MB trigger. 
The RHIC experiments often apply an online timing cut to select events within a 
given primary vertex interval (~ |-z ver tex| <30 cm). Further coincidence with fast 
detectors near mid-rapidity are often also used to suppress background events 
such as beam-gas collisions. Experiment specific event selection is described in 
detail in section 13.21 Figure [9] displays the pseudorapidity coverage of the suite 
of subsystems used to define centrality (both offline and at the trigger level) 
([HTlESlEHlllQD. With the exception of the STAR TPC and FTPCs, all other 
subsystems are intrinsically fast and available for event triggering. 

3.1.2 Centrality Observables In minimum bias p+p and p+p collisions 
at high energy, the charged particle multiplicity dN ev t/dN c h has been measured 
over a wide range of rapidity and is well described by a negative binomial distribu- 
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Figure 10: Charged particle multiplicity (PHOBOS) in Au+Au collisions at 
various center-of-mass energies for < 5 (|57U58p . The different colors represent 
different centrality selections. 



tion (|4ip . However, the multiplicity is also known to scale with the hardness (q 2 ) 
of the collision - the multiplicity for hard jet events is significantly higher than 
that of MB collisions. In heavy ion collisions, we manipulate the fact that the ma- 
jority of the initial state nucleon-nucleon collisions will be analogous to minimum 
bias p+p collisions, with a small perturbation from much rarer hard interactions. 
The final integrated multiplicity of heavy ion events is then roughly described as a 
superposition of many negative binomial distributions, which quickly approaches 
the Gaussian limit. 

N c h can be measured offline by counting charged tracks (e.g, STAR TPC) 
or estimated online from the total energy deposited in a detector divided by 
the typical energy deposition per charged particle (e.g., PHOBOS Paddles). As 
shown in Figured], there is a wide variation in the acceptance of various centrality 
detectors at RHIC. PHOBOS, with the largest acceptance in rj, is well suited to 
measure ^Vi^ (Figure flOj) . These data illustrate two key features of particle 
production in nucleus-nucleus collisions. At a fixed beam energy, there is no 
dramatic change in shape as the centrality changes. However reducing the beam 
energy does change the shape substantially, since the maximum rapidity varies 
as ln(y'sNN)- Thus, the same trigger detector may have a very different overall 
efficiency at different beam energies. 

N c h can be simulated via various methods, but all require the coupling of a 
Glauber calculation to a model of charged particle production, either dynamic 
(e.g. HIJING (|19p ) or static (randomly sampled from a Gaussian, Poisson, or 
negative binomial). Most follow the general prescription that the multiplicity 
scales approximately with A part . For an optical Glauber calculation, simulated 
multiplicity (A^ 111 ) can be calculated semi-analytically assuming that each par- 
ticipant contributes a given value of which is typically drawn from one of the 
aforementioned static probability distributions (42). The same can be done for 
a Monte Carlo Glauber simulation with the added advantage that the detector 
response to such "events" can be simulated, thus enabling an apples-to-apples 
comparison of simulated and measured A c h distributions. Various dynamical 
models of heavy ion collisions exist and can also be coupled to detector simula- 
tions. In all cases, the exact value of A part , N co n, b, and are known for each 
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Figure 11: Spectator energy deposition in the ZDCs as a function of charged 
particle multiplicity in the BBCs (PHENIX). 



event. 

3.1.3 Dividing by Percentile of Total Inelastic Cross Section With 
a measured and simulated d-/V evt /diVch distribution in hand, one can then per- 
form the mapping procedure to extract mean values. Suppose that the measured 
and simulated distributions are both one dimensional histograms. For each his- 
togram, the total integral is calculated and centrality classes are defined in terms 
of fraction of the total integral. Typically the integration is performed from large 
values of to small. For example, the 10-20% most central class is defined by 
boundaries nio and n2o which satisfy 



rmo djV ov t j at- f "20 diV cv t j at 

Joo dAU UJV ch „ „ , Jco dAU UiV< 



r0dA^ dAr -"-land ^ -0.2. (14) 

Joe dN ch UJVch Joo dAT ch UJV ch 

See, for example, Figure The same procedure is performed on both the mea- 
sured and simulated distribution. We note explicitly that n ? ieasured need not equal 
n simuiatod_ T^jg non-trivial fact implies that the mapping procedure is robust to 
an overall scaling of the simulated distribution compared to the measured 
distribution. One can extend the centrality classification beyond one dimension 
by studying the correlation of beam rapidity spectator multiplicity with mid- 
rapidity particle production (Figure [Tl]) . Although the distribution is somewhat 
asymmetric, the naive expectations of Section 13.11 are clearly upheld and the 
mapping procedure proceeds as in the 1-d case. 

Once a centrality class is defined in simulation, the mean values of quantities 
such as iVpart can be calculated for events that fall in that centrality bin. Sys- 
tematic uncertainty in the total measured cross section propagates into a leading 
systematic uncertainty on the Glauber quantities extracted via the mapping pro- 
cess. This uncertainty can be directly propagated by varying the value of the 



Glauber Modeling in Nuclear Collisions 



18 



1 

BBC 




0.2 0.4 0.6 0.8 1.0 

Normalized Si+Tile Multiplicity 



Figure 12: Normalized multiplicity distribution for the Au + 197 Au reaction 
at -y/sNN = 130 GeV measured in the BRAHMS multiplicity array(]44~]). The 
insert shows the correlation pattern of multiplicities measured in the Beam-Beam 
Counter arrays and the multiplicity array. The vertical lines indicate multiplicity 
values associated with the indicated centrality cuts. 

denominator in Equation [14] accordingly and recalculating {N pait ), etc. For ex- 
ample, for Au+Au collisions at ^/iii = 200 GeV in the 10-20% (60-80%) most 
central bin, the STAR collaboration quotes values of (N par t) ~ 234 (21) with 
an uncertainty of ~ 6 (5) from the 5% uncertainty on the total cross section 
alone (|43p . Clearly the uncertainty in total cross section becomes increasingly 
important as one approaches the most peripheral collisions. 

3.2 Experimental Details 

3.2.1 BRAHMS The BRAHMS experiment uses the charged particle mul- 
tiplicity observed in a pseudorapidity range of —2.2 < r/ < 2.2 to determine 
reaction centrality (44,45,46). The multiplicities are measured in a "multiplicity 
array" consisting of an inner barrel of Si strip detectors (SiMA) and an outer 
barrel of plastic scintillator "tile" detectors (TMA) (|37|, l47 p for collisions within 
36 cm of the nominal vertex location. Both arrays cover the same pseudora- 
pidity range for collisions at the nominal vertex. "Beam-Beam" -counter arrays 
(BBC)located on either side of the nominal interaction point at a distance of 
2.2 m of the nominal vertex extend the pseudorapidity coverage for measuring 
charge-particle pseudorapidity densities. These arrays consist of Cherenkov UV- 
transmitting plastic radiators coupled to photomultiplier tubes. Figure [12] shows 
the normalized multiplicity distribution measured in the multiplicity array for 
the 197 Au + 197 Au reaction at y/swN = 130 GeV. The insert shows a correlation 
plot of the multiplicity measured in the Beam-Beam counter array and that in 
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the multiplicity array. The vertical lines indicate the multiplicities corresponding 
to the indicated centrality values. 

The BRAHMS reference multiplicity distribution requires coincident signals 
in the experiment's two ZDC detectors, an interaction vertex located within 30 
cm of the nominal vertex location, and that there be at least four "hits" in 
the TMA. This additional requirement largely removes background contributions 
from beam-residual gas interactions and from very peripheral collisions involving 
only electromagnetic processes. The collision vertex can be determined by either 
the BBC-arrays, the ZDC counters, or a time-projection chamber that is part of 
a mid-rapidity spectrometer arm (TPM1 in Fig. I12|) . 

A simulation of the experimental response based on realistic GEANT3 simu- 
lations (|48p and using the HIJING Monte Carlo event generator (19) for input 
was used to estimate the fraction of the inelastic scattering yield that was missed 
in the experiment's minimum-bias event selection. Multiplicity spectra using the 
simulated events are compared to the experimental spectra. The shapes of the 
spectra are found to agree very well for an extended range of multiplicities in 
the TMA array above the threshold multiplicity set for the event selection. The 
simulated events are then used to extrapolate the experimental spectrum be- 
low the threshold. Using this procedure, it is estimated that the minimum-bias 
event selection criteria selects (93 ± 3)% of the total nuclear cross section for 
Au+Au collisions at ^/snn = 200 GeV, down to (87 ±7)% for Cu+Cu collisions 
at v /snn = 62.4 GeV. 

3.2.2 PHENIX We consider two examples: Au+Au at y/sjm = 200 GeV 
and Cu+Cu at y / SNN = 22.4 GeV. The minimum bias trigger condition for 
Au+Au collisions at -^/snn = 200 GeV was based on Beam-Beam-Counters 
(BBCs) ([38D- The two BBCs (3.1 < \rj\ < 3.9) each consist of 64 photomul- 
tipliers which detect Cherenkov light produced by charged particles traversing 
quartz radiators. On the hardware level a minimum bias event was required to 
have at least 2 photomultiplier hits in each BBC. Some analyses only used events 
with an additional hardware coincidence of the two ZDCs. Moreover, the interac- 
tion vertex along the beam axis {z axis) reconstructed based on the arrival time 
difference in the two BBCs was required to lie within ±30 cm of the nominal 
vertex. 

The efficiency of accepting inelastic Au+Au collisions under the condition of 
having at least two photomultiplier hits in each BBC (Ap^x — 2) was determined 
with the aid of HIJING Monte Carlo events (I19p and a detailed simulation of the 
BBC response (I48p . With an offline vertex cut of ±30 cm these simulations yield 
an efficiency of (92.3±2) % for Au+Au at = 200 GeV (09$. The systematic 

uncertainties reflect uncertainties of (a) dA c h/dy in HIJING, (b) the shape of the 
z-vertex distribution, and (c) the stability of the photomultiplier gains. 

The additional requirement of a ZDC coincidence removes remaining back- 
ground from beam-gas interactions, but possibly also leads to a small inefficiency 
for peripheral collisions. The efficiency of accepting real Au+Au collisions with 
-^pmt — 2 under the condition of a coincidence of the ZDCs was estimated to 
be 99l|'5 Combining the efficiencies of the BBC and ZDC requirement as 
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Figure 13: Measured and simulated distribution of the Pad Chamber 1 hit 
multiplicity used as centrality variable in Cu+Cu collisions at ^/snn = 22.4 GeV 
(PHENIX). 



well as the offline vertex cut PHENIX finds that its sample of minimum bias 
events in Au+Au at ^/snn = 200 GeV corresponds to 91. 4+30 % of the total 
inelastic cross section. Centrality classes were then defined by cuts on the two- 
dimensional distribution of the ZDC energy as a function of the BBC signal as 
shown in Figure [TT1 

As a second example we consider the centrality selection in Cu+Cu collisions 
at ^/snn = 22.4 GeV. At this energy the beam rapidity ?/beam ~ 3.2 lies within 
the pseudorapidity range of the BBCs. The BBCs were still used as minimum 
bias trigger detectors (-ZVpMT > 1 in both BBCs). However, a monotonic relation 
between the BBC signal and the impact parameter was no longer obvious. Thus, 
the hit multiplicity (iVpci) measured with a Pad Chamber detector (150p at mid- 
rapidity (I77I < 0.35) was used as centrality variable. 

The PCI multiplicity distribution was simulated based on a convolution of the 
Apart distribution from Glauber MC and a negative binomial distribution (NBD). 
A non-linear scaling of the average particle multiplicity with A par t was allowed: 
it was assumed that the number of independently decaying precursor particles 
('ancestors', A ances tor) is given by A ancest or = Ng aTt . The number of measured 
PCI hits per precursor particle was assumed follow a NBD: 



P, 



. n 



T(n + k) 



(M/fc) r 



r(n + i)r(fc) (fi/k + i) n + k ' 



(15) 



In a Glauber MC event the NBD was sampled ^ancestor times to obtain the 
simulated PCI multiplicity for this event. The PCI multiplicity distribution 
was simulated for a grid of values for /j,, k , and a in order to find optimal 
values. Figure [13] shows the measured and simulated PCI distribution along 
with the best estimate of the BBC trigger efficiency (cbbc ~ 0.8) corresponding 
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Figure 14: PHOBOS data showing the number of paddle slats hit by charged 
particles, out of 32 total, in triggered events. The data is compared with a full 
Monte Carlo simulation of unbiased HIJING events. The data is normalized 
around N (slats) ~ 16 in order to compare data and MC, to estimate the fraction 
of events lost to the trigger conditions. 

to the difference at small iVpci (see f)49|) for a similar study in Au+Au collisions 
at y / SNN = 19.6 GeV). Given the good agreement between the measured and 
simulated distribution in Figure [13] centrality classes for Cu+Cu collisions at 
snn = 22.4 GeV were defined by identical cuts on the measured and simulated 
iVpci- 

3.2.3 PHOBOS As discussed above, PHOBOS measures centrality with 2 
sets of 16 scintillator paddle counters covering 3 < \n\ < 4.5 ([39]) . Good events are 
defined by having less than 4 ns time difference between the first hit impinging on 
each paddle counter (limiting the vertex range) and either a coincidence between 
the PHOBOS ZDCs or a high summed energy signal in the paddles (to avoid the 
slight inefficiency in the ZDCs at small impact parameter at low energies). 

PHOBOS estimates the observed fraction of the total cross section by mea- 
suring the distribution of the total number of paddle slats, as shown in Fig- 
ure [TH Most of the variation in this quantity essentially measures the very low- 
multiplicity part of the multiplicity distribution, since the bulk of the events have 
sufficient multiplicity to fire all of the paddles simultaneously. The inefficiency is 
determined by matching the "plateau" structure in the data to that in HIJING, 
and measuring how many of the low multiplicity events are missed in the data 
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relative to the MC calculation. This accounts for a variety of instrumental effects 
in an aggregate way, and the difference between the estimated value and 100% 
sets the scale for the systematic uncertainty. 

The total detection efficiency for 200 GeV Au+Au collisions is found to be 
97%, and 88% when requiring more than 2 slats hit on each set of 16 paddle 
counters. This last requirement dramatically reduces background events taken 
to tape, and the relative efficiency is straightforward to measure with the events 
triggered on a coincidence of 1 or more hits in each set of counters. 

For events with much lower multiplicities and/or lower energies, both the mul- 
tiplicity and rapidity reach are substantially smaller. This strongly impacts the 
efficiency of the paddle counters, and thus potentially distorts the centrality esti- 
mation. For these, PHOBOS uses the full distribution of multiplicities measured 
in several ij regions of the Octagon and Ring multiplicity counters, and matches 
them to the distributions measured in a MC simulation, typically HIJING. Once 
the overall multiplicity scale is fit, the difference in the integrals between data 
and MC gives a reasonable estimate of the fraction of observed total cross section. 
3.2.4 STAR STAR defines centrality classes for Cu+Cu and Au+Au (d+Au) 
using charged particle tracks reconstructed in the TPC (FTPC) over full azimuth 
and 1 7/| < 0.5 (2.5 < rj < 4). Background events are removed by requiring the re- 
construction of a primary vertex in addition to either a coincident ZDC (130/200 
GeV Au+Au/Cu+Cu) or BBC (62.4 GeV Au+Au/Cu+Cu) signal. Vertex re- 
construction inefficiency in low multiplicity events reduces the fraction of the 
total measured cross section to, e.g., (95 ± 5)% for 130 GeV Au+Au. For MB 
events, centrality is defined offline by binning the measured dJV cvt /dN^ distribu- 
tion by fraction of total cross section. Glauber calculations are performed using 
a Monte Carlo calculation. STAR enhances central events via an online trigger 
using a coincidence between the MB ZDC condition and large energy deposition 
in the Central Trigger Barrel (CTB), a set of 240 scintillating slats covering full 
azimuth and — 1 < n < 1. After offline cuts the central trigger corresponds to 
the 0-5% most central class of events. STAR has several methods of extracting 
mean values of Glauber quantities. 

(1) STAR reports little dependence on the mean values of N pa _ Tt and N co \\ ex- 
tracted via the aforementioned mapping procedure when vastly different models 
of particle production were used to simulate the charged particle multiplicity. 
Thus, for many analyses (62.4/130/200 GeV Au+Au) STAR bypasses simulation 
of the multiplicity distribution and instead defines centrality bins from the Monte 
Carlo calculated da/dN pa _ rt and dcr/dA^ co ii distributions themselves (H31[5"T |l5"2"l) . 
Mean values of Glauber quantities are extracted by binning the calculated dis- 
tribution (e.g., dcT/dA7p art ) analogously to the measured dAe V t/diV c h. Potential 
biases due to lack of fluctuations in simulated particle production were evaluated 
and found to be negligible for all but the most peripheral events. Further uncer- 
tainties in the extraction of (N pa _ Tt ) and (N co \\) are detailed in reference (f4*3ll5"T1) 
and are dominated by uncertainty in <r^, the Woods Saxon parameters of Au 
and Cu, and the 5% uncertainty in the measured cross section. 

(2) For Cu+Cu and recent studies of elliptic flow fluctuations in Au+Au, STAR 



Glauber Modeling in Nuclear Collisions 



23 



has invoked a full simulation of the TPC multiplicity distribution (|53p . analogous 
to that performed for previous d+Au studies described below. 

(3) For d+Au events, centrality was defined by both measuring and simulating 
the charged particle multiplicity in the FTPC in the direction of the initial Au 
beam. The simulated distribution was constructed using a Monte Carlo Glauber 
model coupled to a random sampling of a NBD. The NBD parameters were 
taken directly from measurements of UA5 collaboration at the same rapidity and 
energy (I41h . For each Monte Carlo event, the NBD was randomly sampled N part 
times. After accounting for tracking efficiency, the simulated distribution was 
found to be in good agreement with the data (|54p . The mean values of various 
Glauber quantities were then extracted as described in Section 13.1.31 A second 
class of events was also used, where a single neutron was tagged in the ZDC in 
the direction of the initial d beam. These "single-neutron" events are essentially 
peripheral p+Au collisions, and the corresponding FTPC multiplicity is again 
well described by the Monte Carlo Glauber based simulation (|54"|) . 

3.3 Acceptance Biases 

Since centrality is estimated using quantities that vary monotonically with par- 
ticle multiplicity, one must be careful to avoid associating fluctuations of an 
observable with fluctuations in the geometric quantities themselves. This is espe- 
cially true when estimating the yield per participant pair, when one is estimating 
the participants from the yield itself. Of course, in heavy ion collisions, the ex- 
traordinarily high multiplicities reduce the effect of auto-correlation bias, as was 
estimated by STAR (|52j). However, the RHIC experiments have found that lower 
multiplicities and lower energies are quite challenging. Estimating the number of 
participants in d+Au proved particularly delicate, due to auto-correlations, which 
were reduced (in HIJING simulations) by using large regions in pseudorapidity 
positioned far forward and backward of mid-rapidity (|7|). 

3.4 Estimating Geometric Quantities 

3.4.1 Total Geometric Cross Section The total geometric cross sec- 
tion for the collision of two nuclei A and B, i.e., the integral of the distributions 
der / db shown in Figure [T5"| is a basic quantity which can be easily calculated in 
the Glauber Monte Carlo approach. It corresponds to all Glauber Monte Carlo 
events with at least one inelastic nucleon-nucleon collision. In ultra-relativistic 
nucleus-nucleus collisions the de Broglie wave length of the nucleons is small com- 
pared to their transverse extent so that quantum mechanical effects are negligible. 
Hence, the total geometric cross section is expected to be a good approximation 
of the total inelastic cross section. For the reaction systems in Figure [T5l (d+Au. 
Cu+Cu, and Au+Au at ^/snn = 200 GeV) the Monte Carlo calculations yield 
ct^+ Au 2180 mb, cj^ +Cu ps 3420 mb, and cr^+ Au « 6840 mb. The systematic 
uncertainties are on the order of 10 % and are dominated by the uncertainties of 
the nuclear density profile. 

Also shown in Figure [15] is a comparison with anf optical calculation of dcr/db, 
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Figure 15: Total geometrical cross section from Glauber Monte Carlo calculations 
(d+Au, Cu+Cu, and Au+Au at - v /Inn = 200 GeV). The dashed line represents 
an optical limit calculation for Au+Au. 



which shows the effect described in Section 12,31 Optical limit calculations do 
not naturally contain the terms in the multiple-scattering integral where nu- 
cleons "hide" behind each other. This leads to a slightly larger cross section 
( a geo^pticai ~ 7280 mb). While this seems like a small perturbation on a^ B , it 
has a surprisingly large effect on the extraction of N paxt and iV co n . This does not 
come from any fundamental change in the mapping of impact parameter onto 
those variables. Figure [5] shows the mean value of N co \i (upper curve) and -/V part 
(lower curve) as a function of b, where it is seen that the two track each other very 
precisely over a large range in impact parameter, well within the range usually 
measured by the RHIC experiments. The problem comes in when dividing a sam- 
ple up into bins in fractional cross section. While it is straightforward to estimate 
the most central bins, one finds a systematic difference of -/Vp art between the two 
calculations as the geometry gets more peripheral. This is shown in Figure [16] 
for y / SNN = 130 GeV by comparing the Monte Carlo calculation in HIJING with 
the optical limit calculation in (|42|) . 

3.4.2 Participants (N pait ) and Binary Collisions (-/V co ii) As described 
in Section [21 the Glauber model is a multiple collision model which treats a 
nucleus-nucleus (A+B) collision as an independent sequence of nucleon-nucleon 
collisions. A participating nucleon or wounded nucleon is defined as a nucleon 
which undergoes at least one inelastic nucleon-nucleon collision. The centrality 
of a A+B collision can be characterized both by the number of participating 
nucleons (-/Vp art ) and by the number of binary nucleon-nucleon collisions (N co u). 
The average number of participants (N paTt ) and nucleon-nucleon collisions {N co \\) 
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Figure 16: The ratio of N pai t calculated with a Glauber Monte Carlo (N^{) to 
that calculated with an optical approximation (-/Vp^ cal ), for the same fraction of 
the total inelastic Au+Au cross section (■ v /snn = 130 GeV), plotted as a function 



as a function of the the impact parameter b are shown in Figure [T7] for Au+Au 
and Cu+Cu collisions at -^snn = 200 GeV. The event-by-event fluctuations of 
these quantities for a fixed impact parameter are illustrated by the scatter plots. 

The shapes of the iV part and N co \\ distributions shown in Figure [TH] for Au+Au 
collisions reflect the fact that peripheral nucleus-nucleus collisions are more likely 
than central collisions. {N pSLTt ) and {N co \\) for a given experimental centrality 
class, e.g., the 10% most central collisions, depend on the fluctuations of the 
centrality variable which is closely related to the geometrical acceptance of the 
respective detector. By simulating the fluctuations of the experimental centrality 
variable and applying similar centrality cuts as in the analysis of real data one ob- 
tains iVpart and N co \\ distribution for each centrality class. For peripheral classes 
the bias introduced by the inefficiency of the experimental minimum bias trigger 
needs to be taken into account by applying a corresponding trigger threshold on 
the Glauber Monte Carlo events. Experimental observables like particle multi- 
plicities can then be plotted as a function of the mean value of N pil n and N co \\ 
distributions. 

The systematic uncertainties of iV part , N co u, and other calculated quantities are 
estimated by varying various model parameters. Figure [19] shows such a study 
for Au+Au collisions at ^/snn = 200 GeV (PHENIX). The following effects were 
considered: 

• The default value of the nucleon-nucleon cross section of <tnn = 42 mb was 
changed to 39 mb and 45 mb. 

• Woods-Saxon parameters were varied to determine uncertainties related to 
the nuclear density profile. 

• Effects of a nucleon hard core were studied by requiring a minimum distance 
of 0.8 fm between two nucleons of the same nucleus without distorting the 
radial density profile. 

• Parameters of the BBC and ZDC simulation (e.g. parameters describing 
the finite resolution of these detectors) were varied. 
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Figure 17: Average number of participants ((iVp ar t)) and binary nucleon-nucleon 
collisions ((iV co ii)) along with event-by-event fluctuation of these quantities in the 
Glauber Monte Carlo calculation as a function of the impact parameter b. 




Figure 18: iVp ar t and N co \\ distribution in Au+Au collisions at y / SNN = 200 GeV 
from a Glauber Monte Carlo calculation. By applying cuts on simulated centrality 
variables, in this case the BBC and ZDC signal as measured by PHENIX, one 
obtains iVp art and iV co ii distributions for the respective centrality class. 
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Figure 19: Effect of various parameters in the Glauber Monte Carlo calculation 
on iVp art and A^jj for Au+Au collisions at ^/snn = 200 GeV. 
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• The black disk nucleon-nucleon overlap function was replaced by "gray disk" 
and Gaussian overlap function (13T1) without changing the total inelastic 
nucleon-nucleon cross-section. 

• The origin of the centrality cuts applied in the scatter plot of ZDC vs. BBC 
space was modified in the Glauber calculation. 

• The uncertainty of the efficiency of the minimum bias trigger leads to uncer- 
tainties as to which percentile of the total inelastic cross section actually is 
selected with certain centrality cuts. The centrality cuts applied on the cen- 
trality observable simulated with the Glauber MC were varied accordingly 
to study the influence on (N pa _ Tt ) and (iV co ii). 

• Even if the minimum bias trigger efficiency were precisely known potential 
instabilities of the centrality detectors could lead to uncertainties as to 
which percentile of the total cross section is selected. This has been studied 
by comparing the number of events in each experimental centrality class for 
different run periods. The effect on (iV par t) and (iV co u) was again studied 
by varying the cuts on the simulated centrality variable accordingly. 

The total systematic uncertainty indicated by the shaded boxes in Figure [19] 
were obtained by adding the deviations from the default result for each of the 
items in the above list in quadrature. The uncertainty of N pavt decreases from 
~ 20 % in peripheral collisions to ~ 3 % in central Au+Au collisions. N co n has 
similar uncertainties as iVp ar t for peripheral Au+Au collisions. For N part > 100 
(or N co \\ > 200) the relative systematic uncertainty of N co u remains constant at 
about ~ 10%. Similar estimates for the systematic uncertainties of N paTt and 
TVcon at the CERN SPS energy of v /snn = 17.2 GeV were reported in ([5%]). 

For the comparison of observables related to hard processes in A+A and p+p 
collisions it is advantageous to introduce the nuclear overlap function {Ta,b}{ for 
a certain centrality class f (see section 14. 2\i which is calculated in the Glauber 
Monte Carlo approach as 



The uncertainty of the inelastic nucleon-nucleon cross section cr^2 doesn't con- 
tribute to the systematic uncertainty. Apart from this (Tab){ has the same 
systematics uncertainties as (N co u)f. 

3.4.3 Eccentricity One of the surprising features of the RHIC data was 
the strong event-by-event asymmetries in the azimuthal distributions. This has 
been attributed to the phenomenon of "elliptic flow" , the transformation of spa- 
tial asymmetries into momentum asymmetries by means of hydrodynamic evolu- 
tion. For any hydrodynamic model to be appropriate, the system must be suf- 
ficiently opaque (where opacity is the product of density times interaction cross 
section) such that the system equilibrates locally at early times. This suggests 
that the relevant geometric quantity for controlling elliptic flow is the "shape" of 
the overlap region, which sets the scale for the gradients that develop. 

The typical variable used to quantify this shape is the "eccentricity", defined 




(16) 



Glauber Modeling in Nuclear Collisions 



29 




Figure 20: Calculation of eccentricity in Au+Au and Cu+Cu collisions as a 
function of -/Vp art . Both standard and participant eccentricity are shown. 



as 

(Y 2 -X 2 ) , . 

Just as with other variables discussed here, there are two ways to calculate this. 
In the optical limit, one performs the averages at a fixed impact parameter, 
weighting by either the local participant or binary collision density. In the Monte 
Carlo approach, one simply calculates the moments of the participants them- 
selves. Furthermore, one can calculate these moments with the X axis oriented 
in two natural frames. The first is along the nominal reaction plane (estimated 
using spectator nucleons). The second is along the short principal axis of the 
participant distribution itself (1561) . The only mathematical difference between 
the two definitions involves the incorporation of the correlation coefficient {XY): 



fcpart — 5 5 (18) 

A comparison of the two definitions is shown in Figure [20j One sees very 
different limiting behavior at very large and small impact parameter. At large 
impact parameter, fluctuations due to small numbers of participants drive e s td — ► 
0, but e part — ► 1. As b = 0, e st d also goes to zero as the system becomes radially 
symmetric while e pa rt now picks up the fluctuations and remains finite. The 
relevance of these two quantities to actual data will be discussed in sections 14.31 
andEOJ 



Glauber Modeling in Nuclear Collisions 



30 




Figure 21: Total inclusive charged-particle multiplicity (N,±) divided by iV par t/2 
from PHOBOS data at four RHIC energies. The data is compared with p + p 
data or interpolations to unmeasured energies at N paTt = 2. 



4 Geometric aspects of p+A and A+A phenomena 

4.1 Inclusive Charged-Particle Yields (total and mid-rapidity) 

The total multiplicity in hadronic reactions is a measure of the degrees of free- 
dom released in the collision. In the 1970s it was found that the total number 
of particles produced in proton-nucleus (p + A) collisions was proportional to 
the number of participants, i.e. oc N pa _ Tt = v + 1, where v is defined as the 
number of struck nucleons in the nucleus (16ip . This experimental fact was instru- 
mental in establishing -/V part as a fundamental physical variable. The situation 
became more interesting when the total multiplicity was measured in Au+Au at 
4 RHIC energies, spanning an order of magnitude in y / SNN, and was found to be 
approximately proportional to N paTt there as well. This is shown in Fig]21l with 
PHOBOS data from Refs. (|58| 59) and is striking if one considers the variety of 
physics processes that should contribute to the total multiplicity. 

By contrast, the inclusive charged-particle density near mid-rapidity (dN / dr)(\rf\ < 
1)) does not scale linearly with N paT t/2. This is shown in Fig. [22] with STAR data, 
with ./Vp ar t estimated both from an optical calculation (left) as well as a Glauber 
Monte Carlo (right) (|62p . The comparison shows why care must be taken in the 
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Figure 22: Inclusive charged particle multiplicity near midrapidity (dN/drj(\rj\ < 
1)) divided by the number of participating nucleon pairs (A^ part /2) estimated 
using the optical approximation (left) and a Glauber Monte Carlo (right) from 
STAR. The data is compared with two-component fits and a calculation based 
on parton saturation. 



estimation of -/Vp ar t, since using one or the other gives better agreement with very 
different models. The saturation model of Eskola et al (f63|) does show scaling 
with iVp ar t and agrees with the data if -/V part is estimated with an optical cal- 
culation. However, it disagrees with the data when the GMC approach is used. 
Better agreement with the data can be found with a so-called "two component" 
model, e.g. Ref. (02]): 

— = npp [(1 - x)N paTt + xN coU ] (19) 

This model can be fit to the data by a judicious choice of x, the parameter 
which controls the admixture of "hard" particle production. However there is no 
evidence of any energy dependence to this parameter (65,64) from y / s/viv = 19.6 
to 200 GeV, suggesting that the source of this dependence has little to nothing 
to do with hard or semi-hard processes at all. 



4.2 Hard Scattering: T AB Scaling 

The number of hard processes between point-like constituents of the nucleons in 
a nucleus-nucleus collision is proportional to the nuclear overlap function Tab(&) 
([6^ [67 1 [68 | [69 l [70 l[7T]) . This follows directly from the factorization theorem in 
the theoretical description of hard interactions within perturbative QCD (|72h . 
In detail, the average yield for a hard process with cross section cr?;L d in p+p 
collisions per encounter of two nuclei A and B with impact parameter b is given 
by 
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Here Tab(^) is normalized so that J TAB(b) d 2 b = AB. c^ard can > e -5-' represent 
the cross-section for the production of charm-anticharm quark-pairs (cc) or high- 
Pt direct photons in proton-proton collisions. 

At large impact parameter not every encounter of two nuclei A and B leads to 
an inelastic collision. Hence, the average number of hard processes per inelastic 
A+B collision is given by 



= ^1 -hard (21) 



where P\^^(b) is the probability of an inelastic A+B collision. In the optical 
limit pf^^(b) is given by (cf. Eq. [5]) 

A+B / t\ -i ( -i ^NN^AB(^)\ AB 
Pinel W = 1 - ( 1 - ^incl AR J ■ (22) 

Particle yields at RHIC are usually measured as a function of the transverse 
momentum (pt)- If an invariant cross section dcr pp /dpr for a hard scattering 
process which leads to the production of a certain particle x has been measured 
in p+p collisions, then according to Eq. [21] the invariant multiplicity of x per 
inelastic A+B collisions with impact parameter b is given by 

1 d7V A + B _ T A b (b) da™ 



This baseline expectation is purely based on nuclear geometry and assumes the 
absence of any nuclear effects. In reality one needs to average over a certain 
impact parameter distribution. As an example we consider a centrality class f 
which corresponds to a fixed impact parameter range 61 < b < b<i- Taking the 
average in this range over the impact parameter distribution (weighting factor 
da AB /d6 = 2tt&p a + b (&)) leads to 



1 diV A+B 



iV AB d PT 

with 



(T AB )f • ^2- (24) 



, M / db 2tt6 Tab b 

/ f d 2 bT AB (b) _ by ( s 

( AB)f = fd%v A+B (b) ~ ^ ' ( ] 

Jfd b Pind W /d6 2vr6p A + B (6) 
61 

Averaging over the full impact parameter range (61 = 0, 62 = 00) yields (TX.B)f = 
AB / a^ B . In the Glauber Monte Carlo approach (T\B)f for a certain centrality 
class is calculated as 

(T AB )f = (JWf/otad (26) 

where the averaging is done for all A+B collisions with at least one inelastic 
nucleon-nucleon collision and whose simulated centrality variable belongs to cen- 
trality class f. 
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Figure 23: Nuclear modification factor Raa in Au+Au collisions at y / SNN = 
200 GeV for direct-photon and neutral-pion yields integrated above = 
6 GeV/c. The dashed lines indicate the systematic uncertainties of (TAB)f used 
in the calculation of Raa- 

In order to quantify nuclear effects on particles production in hard scattering 
processes the nuclear modification factor -Rab(pt) is defined as 



At high (pt > 2 — 3 GeV/c for hadrons and pp > 4 — 6 GeV/c for direct 
photons) particle production is expected to be dominated by hard processes such 
that, in the absence of nuclear effects, Rab should be unity. Due to their elec- 
tromagnetic nature high-px direct photons are essentially unaffected by the hot 
and dense medium produced in a nucleus-nucleus collisions so that they should 
exhibit Tab scaling. Fig. [23] shows the nuclear modification factor for direct- 
photon and neutral-pion yields integrated above pt = 6 GeV/c as a function of 
A^art (73). Direct photons indeed follow Tab scaling over the entire centrality 
range whereas neutral pions are strongly suppressed in central collisions. This 
is one of the major discoveries at RHIC. The direct-photon measurement is an 
experimental proof of Tab scaling of hard processes in nucleus-nucleus collisions. 
With this observation the most natural explanation for the suppression of high- 
Pt neutral pions is energy loss of partons from hard scattering in a quark-gluon 
plasma (jet-quenching) f^ElEKH]) . 

4.3 Eccentricity and relation to Elliptic Flow 

Hydrodynamic calculations suggest that spatial asymmetries in the initial state 
are mapped directly into asymmetries in the final state momentum distribu- 
tion. At mid-rapidity, these asymmetries are manifest in the azimuthal (</>) dis- 
tributions of inclusive and identified charged particles, with the modulation of 




(27) 
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dN/dcfr ~ 1 + 2t>2Cos[2(0 — ^r)] characterized by the 2 nd Fourier coefficient v?, 
where ^r defines the angle of the reaction plane for a given event. It is generally 
assumed that t>2 is proportional to the event eccentricity e, which was introduced 
previously. Glauber calculations are used to estimate the eccentricity, either for 
an ensemble of events or on an event-by-event basis. 

For much of the RHIC program, both calculations were typically carried out in 
the "standard" reference frame, with the X-axis oriented along the reaction plane. 
Using this calculation method was apparently sufficient to compare hydrodynamic 
calculations with Au+Au data. However, it was always noticed that the most 
central events, which should trend to e = tended to have a significant V2 value. 
This led to the study of the "participant eccentricity" , calculated with the X axis 
oriented along the short principal axis of the approximately-elliptical distribution 
of participants in a Monte Carlo approach (|56|) . described above. 

Fig. [241 from Ref. ([741) shows V2, the second Fourier coefficient ((cos(2[^>— 
of the inclusive particle yield relative to the estimated reaction plane angle, as a 
function of -/Vp art . In hydrodynamic models, V2 is proportional to the eccentricity, 
suggesting that V2/e should be a scaling variable. While the raw values of V2 as 
a function of N paiTt peak at similar levels in Au+Au and Cu+Cu, it is found that 
dividing by the standard eccentricity makes the two data sets diverge. However, 
dividing by e par t shows that the two systems have similar i>2/e P art at the same 
Apart- This shows that the participant eccentricity, a quantity calculated in a 
simple Glauber Monte Carlo approach, drives the hydrodynamic evolution of the 
system for very different energy and geometries. 

4.4 Eccentricity Fluctuations 

As described above, one of the most spectacular measurements at RHIC was the 
large value of elliptic flow in Au+Au collisions, suggestive of a "Perfect Liquid." 
After the initial measurement (I75p . much attention was given to potential biases 
that could artificially inflate the extraction of V2 from the data, such as "non-flow" 
effects (e.g., correlations from jet fragmentation, resonance decay) and event-by- 
event fluctuations in V2 itself. Reference (|32|) was one of the first analyses to study 
the effects of fluctuations on extraction of V2- Using the assumption V2 ~ e s td, 
fluctuations in e st d were studied using a Monte Carlo Glauber calculation and 
comparing (e" td ) to (e st d)™. Fluctuations were found to play a significant role, 
where different methods of extraction (e.g. 2-particle vs. higher order cumulants) 
gave results differing by as much as a factor of two, with the most significant 
differences found for the most central (0-5%) and most peripheral (60-80%) events 
classes. 

Recently in references (|53[) and (|76p . the STAR and PHOBOS collaborations 
have reported measurements of not only the (1J2), but also the r.m.s. width 
u v% . Figure [25] shows the distribution of cr V2 /(v2) vs (b) from Au+Au data. The 
measurements are compared to calculations from various dynamical models as 
well as Monte Carlo Glauber calculations using both <T £atd /(e s td) and o"e part /(epart}- 
Clearly the e st d description is ruled out while the e pa rt description is in good 
agreement within the measured uncertainties, implying that the measurements 
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Figure 24: (a) PHOBOS Measurements of v 2 for 62.4 and 200 GeV Cu+Cu and 
Au+Au collisions (four system/energy combinations in all, from Ref. (|74D ). (b) 
V2 divided by standard eccentricity and (c) v 2 divided by participant eccentricity, 
showing an approximate scaling for all the system/energy combinations. 
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Figure 25: (a) Elliptic flow mean value (top) and r.m.s. width scaled by the mean 
(bottom) as measured by STAR, compared to various dynamical models and a 
Monte Carlo Glauber calculation ([53]) . (b) Gaussian width of Vi divided by mean 
v 2 from PHOBOS d7BJ. 

are sensitive to the initial conditions. The agreement with the €p ar t Glauber 
calculation further implies that the measured V2 fluctuations are fully accounted 
for by the fluctuations in the initial geometry, leaving little room for other sources 
(e.g., Color Glass Condensate). We note that these analyses are new and the 
physics conclusions are far from final. However, this is another excellent example 
where Glauber calculations are critical in interpreting RHIC data. 

4.5 J/ip absorption in normal nuclear matter 

Due to the large mass of the charm quark cc pairs are expected to be produced 
only in hard processes in the initial phase of a nucleus-nucleus collision. The pro- 
duction rate for cc pairs is thus calculable within perturbative QCD which makes 
them a calibrated probe of later stages of a heavy ion collision. In particular, it 
was suggested that free color charges in a quark-gluon plasma could prevent the 
formation of a J/tp from the initially produced cc pairs so that J ftp suppression 
was initially considered a key signature of a QGP formation (77). However, a 
suppression of J/ip's relative to the expected production rate for hard processes 
was already seen in proton-nucleus (p+A) collisions. Thus, it became clear that 
the "conventional" J ftp suppression in p+A collisions needed to be quantified 
and extrapolated to A+B collisions before any conclusions could be drawn about 
a possible QGP formation in A+B collisions. Cold nuclear matter effects which 
affect J 7V> production include the modification of the parton distribution in the 
nucleus (shadowing) and the absorption of pre-resonant cc pairs (78). For the 
extrapolation of the latter effect from p+A to A+A collisions the Glauber model 
is frequently used (f79l) . as described in the following. 

We first concentrate on p+A collisions (left panel of Fig. [26]) . Conventional J/ip 
suppression is thought to be related to the path along the z axis in normal nuclear 




Figure 26: Sketch illustrating the calculation of the J/tfj absorption in normal 
nuclear matter in proton-nucleus (left panel) and nucleus-nucleus (right panel) 
collisions). The dashed arrows indicate the paths of a cc pair created in a hard 
process. 



matter a pre-resonant cc pair created at point (b, za) needs to travel. Since the cc 
pair is created in a hard process the location of its production is indeed rather well 
defined. Usually the effects of processes which inhibit the formation of a J ftp from 
the pre-resonant cc pair are parameterized with a constant "absorption" cross 
section <r a b s - Furthermore, formation time effects are often neglected. Under 
these assumptions the probability for the break-up of a pre-resonant cc state 
created at (b, Za) in a collision with a certain nucleon j of nucleus A then reads 

oo 

Pabs(b,^A) = 0"abs^A>(b,<2A) witn T A > (b, Z A ) = ^ /SaO 3 , z) dz . (28) 

za 

Here pa is the density profile of nucleus A normalized so that integration over 
full space yields unity. Hence, the total survival probability for the cc pair is 

p^ rv (b, z A ) = (l - o"abs TA>(b, za)) ~ exp (-(yl - 1) a abs f A> (h, z A )) 

(29) 

where the approximation holds for large nuclei 4>1. The term A — 1 reflects 
the fact that the cc-producing nucleon doesn't contribute to the absorption. The 
spatial distribution of the produced cc pairs follows pA(b,z). Thus, for impact 
parameter averaged p+A collisions ("minimum bias") the expression for the J /ip 
absorption in the Glauber model reads 

S p+ a = = J d2bdz A PA(b, z A ) pi v (b, z A ) . (30) 

Fitting the Glauber model expectation to p+A data yields absorption cross- 
sections on the order of a few mb at CERN SPS energies (80). 

As illustrated by the dashed arrows in the right panel of Fig. 1261 a cc pair 
created in a collision of two nuclei A and B has to pass through both nuclei. 
Analogous to Eq. [28] and [29] the J/ ip survival probability for the path in nucleus 
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B can be written as 

2B 



,(s-b, z B ) = (l - o-abs Tb<(s - b, z B )) , T B <(s-b, z B ) = J p B (s-b, z) dz 



(31) 

The spatial cc production probability density follows pa(s, z\) ■ p B (s — b, z B ) so 
that the normal J/^ suppression as estimated with the Glauber model for A+B 
collisions with fixed impact parameter b is given by 

dS A+B ( , 1 da(AB -+ J/VQ 

d 2 6 l j ABa(p + p^ J/4>) d 2 6 1 J 

= y d 2 s dz A dz B p A (s, z A )p B (s - b, z B ) p^ rv (s, z A )pf nrv (s - b, z B ) 

Here dS < A+B/d 2 6 is normalized so that / d 2 6dS , A+ B (b)/d 2 6 = 1 for <7 a b s = 0. 
Sometimes the J/tp suppression observed in p+A collisions is extrapolated to 
A+B collisions within the Glauber framework by calculating an effective path 
length L = L A + L B so that the expected normal J/ij) suppression can be written 

as cgisaETj 

Sa+b = exp(-L po cr abs ) (33) 

where po = 0.17 fm _1 is the nucleon density in the center of heavy nuclei. At 
the CERN SPS a suppression stronger than expected from the absorption in cold 
nuclear matter has been observed in central Pb+Pb collisions (|81h . This so-called 
"anomalous" J /if) suppression has been discussed as a potential signal for a QGP 
formation at the CERN SPS energy. 
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5 Discussion and the Future 

The Glauber model as used in ultra-relativistic heavy-ion physics is purely based 
on nuclear geometry. What is left from its origin as a quantum mechanical 
multiple scattering theory is the assumption that a nucleus-nucleus collisions can 
be viewed as sequence of nucleon-nucleon collisions and that individual nucleons 
travel on straight line trajectories. With the number of participants -/Vpart and the 
number of binary nucleon-nucleon collisions N co \\ the Glauber model introduces 
quantities which are essentially not measurable. Only the forward energy in 
fixed-target experiments has a rather direct relation to N pa _ Tt . 

The motivation for the use of these rather theoretical quantities is manifold. 
One of the main reasons for the use of geometry-related quantities like N paT t 
calculated with the Glauber model is the possibility to compare centrality depen- 
dent observables measured in different experiments. Moreover, the comparison 
of different reaction system as a function of geometric quantities often leads to 
new insights. Basically all experiments calculate N part and N co \\ in a similar way 
using a Monte Carlo implementation of the Glauber model so that the theoretical 
bias introduced in the comparisons is typically small. Thus, the Glauber model 
provides a crucial interface between theory and experiment. 

The widespread use of the Glauber model is related to the fact that indeed 
many aspects of ultra-relativistic nucleus-nucleus collisions can be understood 
purely based on geometry. A good example is the total charged particle mul- 
tiplicity which scales as iVp art over a wide centrality and center-of-mass energy 
range. Another example is the anisotropic momentum distribution of low-pr 
(pt <~ 2 GeV/c) particles with respect to the reaction plane. This so-called el- 
liptic flow has its origin in the spatial anisotropy of the initial overlap volume in 
non-central nucleus-nucleus collisions. It is a success of the Glauber model that 
event-by-event fluctuations of the spatial anisotropy of the overlap zone as calcu- 
lated in the Monte Carlo approach appear to be relevant for the understanding 
of the measured elliptic flow. In this way, a precise understanding of the Glauber 
picture has been of central concern for understanding the matter produced at 
RHIC as a "near perfect fluid". 

The study of particle production in hard scattering processes is another im- 
portant field of application for the Glauber model. According to the QCD fac- 
torization theorem the only difference between p+p and A+A collisions in the 
perturbative QCD description in the absence of nuclear effects is the increased 
parton flux. This corresponds to a scaling of the invariant particle yields with 
the number of binary nucleon-nucleon collisions (AT coll ) as calculated with the 
Glauber model. The scaling of hard processes with iV co ii or Tab was confirmed 
by the measurement of high-pp direct photons in Au+Au collisions at RHIC. 
This supported the interpretation of a deviation from Tab scaling for neutral pi- 
ons and other hadrons (high-px hadron suppression) as a result of parton energy 
loss in a quark-gluon plasma. 

Future heavy ion experiments, both at RHIC and at the LHC will further push 
our understanding of nuclear geometry. As RHIC experiments study more com- 
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plex multiparticle observables, the understanding of fluctuations and correlations 
even in something as apparently simple as the Glauber Monte Carlo will become 
a limiting factor in interpreting data. And as the study of high px phenomena 
involving light and heavy flavor becomes prominent in the RHIC II era, the un- 
derstanding of nuclear geometry, both experimentally and theoretically, will limit 
the experimental systematic errors. At the LHC, the precision of the geometric 
calculations will be limited by the knowledge of cj^, which should be measured 
in the first several years of the p + p program. After that, Glauber calculations 
will be a central part of understanding the baseline physics of heavy ions at the 
LHC in terms of nuclear geometry. It is hoped that this review will prepare the 
next generation of relativistic heavy ion physicists for tackling these issues. 
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